library(stringr)
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.4 ✓ purrr 0.3.4
✓ tibble 3.1.2 ✓ dplyr 1.0.7
✓ tidyr 1.1.3 ✓ forcats 0.5.1
✓ readr 1.4.0
── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x readr::col_factor() masks scales::col_factor()
x purrr::discard() masks scales::discard()
x tidyr::expand() masks Matrix::expand()
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
x tidyr::pack() masks Matrix::pack()
x tidyr::unpack() masks Matrix::unpack()
response <- try(system('~/google-cloud-sdk/bin/gcloud projects list --quiet', intern = T))
projectid <- strsplit(response[2], " ")[[1]][1]
options(na.action = "na.fail")
source("./dredge_functions.R")
dredge_and_subset <- function(data) {
model <- lm(
presence_ratio ~ foraging_niche + trophic_niche + is_noctural + pc1 + pc2 + pc3 + pc4,
data=data
)
dredge_result <- dredge(model)
summary(model.avg(dredge_result, subset = delta < 10))
}
load_species_data <- function(poolname) {
filename <- str_replace('species_metrics_##POOL_NAME##.csv', '##POOL_NAME##', poolname)
if (!file.exists(filename)) {
column_name = str_replace_all('present_in_##POOL_NAME##_pool', '##POOL_NAME##', poolname)
if (poolname == 'both') {
column_name = 'present_in_both_pools'
}
sql <- str_replace_all("
WITH species AS (
SELECT
scientific_name,
count(*) AS regional_pool_count,
countif(present_in_city = true) AS city_count
FROM model.regional_species_filtered
JOIN model.taxonomy USING (scientific_name)
WHERE ##COL_SELECT## = true
GROUP BY scientific_name
)
SELECT
species.*,
body_morphspace.pc1,
body_morphspace.pc2,
body_morphspace.pc3,
body_morphspace.pc4,
trophic_niche,
foraging_niche,
is_noctural,
ttc.cluster,
ttc.cluster_half,
ttc.cluster_double
FROM species
JOIN model.taxonomy USING (scientific_name)
JOIN datasci.taxonomic_trait_clusters ttc USING (scientific_name)", '##COL_SELECT##', column_name)
tb <- bq_project_query(projectid, sql)
data <- bq_table_download(tb)
write_csv(data, filename)
}
data <- read_csv(filename)
data[is.na(data$foraging_niche),]$foraging_niche <- 'Omnivore'
data$foraging_niche = as.factor(data$foraging_niche)
data$trophic_niche = as.factor(data$trophic_niche)
data$is_noctural = as.factor(data$is_noctural)
data$presence_ratio = data$city_count / data$regional_pool_count
data
}
merlin_data <- load_species_data('merlin')
── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
scientific_name = col_character(),
regional_pool_count = col_double(),
city_count = col_double(),
pc1 = col_double(),
pc2 = col_double(),
pc3 = col_double(),
pc4 = col_double(),
trophic_niche = col_character(),
foraging_niche = col_character(),
is_noctural = col_logical(),
cluster = col_double(),
cluster_half = col_double(),
cluster_double = col_double()
)
merlin_data
merlin_result <- dredge_and_subset(merlin_data)
Fixed term is "(Intercept)"
merlin_summ <- model_summary('merlin', 'species', merlin_result)
merlin_summ
birdlife_data <- load_species_data('birdlife')
── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
scientific_name = col_character(),
regional_pool_count = col_double(),
city_count = col_double(),
pc1 = col_double(),
pc2 = col_double(),
pc3 = col_double(),
pc4 = col_double(),
trophic_niche = col_character(),
foraging_niche = col_character(),
is_noctural = col_logical(),
cluster = col_double(),
cluster_half = col_double(),
cluster_double = col_double()
)
birdlife_result <- dredge_and_subset(birdlife_data)
Fixed term is "(Intercept)"
birdlife_summ <- model_summary('birdlife', 'species', birdlife_result)
birdlife_summ
both_data <- load_species_data('both')
── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
scientific_name = col_character(),
regional_pool_count = col_double(),
city_count = col_double(),
pc1 = col_double(),
pc2 = col_double(),
pc3 = col_double(),
pc4 = col_double(),
trophic_niche = col_character(),
foraging_niche = col_character(),
is_noctural = col_logical(),
cluster = col_double(),
cluster_half = col_double(),
cluster_double = col_double()
)
both_result <- dredge_and_subset(both_data)
Fixed term is "(Intercept)"
both_summ <- model_summary('both', 'species', both_result)
both_summ
either_data <- load_species_data('either')
── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
scientific_name = col_character(),
regional_pool_count = col_double(),
city_count = col_double(),
pc1 = col_double(),
pc2 = col_double(),
pc3 = col_double(),
pc4 = col_double(),
trophic_niche = col_character(),
foraging_niche = col_character(),
is_noctural = col_logical(),
cluster = col_double(),
cluster_half = col_double(),
cluster_double = col_double()
)
either_result <- dredge_and_subset(either_data)
Fixed term is "(Intercept)"
either_summ <- model_summary('either', 'species', either_result)
either_summ
all_species_results <- full_join(full_join(merlin_summ, birdlife_summ), full_join(both_summ, either_summ))
Joining, by = c("explanatory", "model")
Joining, by = c("explanatory", "model")
Joining, by = c("explanatory", "model")
write_csv(all_species_results, "species_result.csv")


library(ggpubr)

jpeg(width = 1024, height = 1024)
myplot
dev.off()
null device
1
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQpgYGB7cn0KbGlicmFyeShzdHJpbmdyKQpsaWJyYXJ5KHRpZHl2ZXJzZSkKYGBgCgpgYGB7cn0KcmVzcG9uc2UgPC0gdHJ5KHN5c3RlbSgnfi9nb29nbGUtY2xvdWQtc2RrL2Jpbi9nY2xvdWQgcHJvamVjdHMgbGlzdCAtLXF1aWV0JywgaW50ZXJuID0gVCkpCnByb2plY3RpZCA8LSBzdHJzcGxpdChyZXNwb25zZVsyXSwgIiAiKVtbMV1dWzFdCmBgYAoKYGBge3J9Cm9wdGlvbnMobmEuYWN0aW9uID0gIm5hLmZhaWwiKSAKYGBgCgpgYGB7cn0Kc291cmNlKCIuL2RyZWRnZV9mdW5jdGlvbnMuUiIpCmBgYAoKYGBge3J9CmRyZWRnZV9hbmRfc3Vic2V0IDwtIGZ1bmN0aW9uKGRhdGEpIHsKICBtb2RlbCA8LSBsbSgKICAgIHByZXNlbmNlX3JhdGlvIH4gZm9yYWdpbmdfbmljaGUgKyB0cm9waGljX25pY2hlICsgaXNfbm9jdHVyYWwgKyBwYzEgKyBwYzIgKyBwYzMgKyBwYzQsCiAgICBkYXRhPWRhdGEKICApCiAgZHJlZGdlX3Jlc3VsdCA8LSBkcmVkZ2UobW9kZWwpCiAgc3VtbWFyeShtb2RlbC5hdmcoZHJlZGdlX3Jlc3VsdCwgc3Vic2V0ID0gZGVsdGEgPCAxMCkpCn0KYGBgCgpgYGB7cn0KbG9hZF9zcGVjaWVzX2RhdGEgPC0gZnVuY3Rpb24ocG9vbG5hbWUpIHsKICBmaWxlbmFtZSA8LSBzdHJfcmVwbGFjZSgnc3BlY2llc19tZXRyaWNzXyMjUE9PTF9OQU1FIyMuY3N2JywgJyMjUE9PTF9OQU1FIyMnLCBwb29sbmFtZSkKICAKICBpZiAoIWZpbGUuZXhpc3RzKGZpbGVuYW1lKSkgewogICAgY29sdW1uX25hbWUgPSBzdHJfcmVwbGFjZV9hbGwoJ3ByZXNlbnRfaW5fIyNQT09MX05BTUUjI19wb29sJywgJyMjUE9PTF9OQU1FIyMnLCBwb29sbmFtZSkKICAgIGlmIChwb29sbmFtZSA9PSAnYm90aCcpIHsKICAgICAgY29sdW1uX25hbWUgPSAncHJlc2VudF9pbl9ib3RoX3Bvb2xzJwogICAgfQogICAgCiAgICBzcWwgPC0gc3RyX3JlcGxhY2VfYWxsKCIKICAgICAgICAgV0lUSCBzcGVjaWVzIEFTICgKICAgICAgICAgICAgICBTRUxFQ1QgCiAgICAgICAgICAgICAgICAgIHNjaWVudGlmaWNfbmFtZSwKICAgICAgICAgICAgICAgICAgY291bnQoKikgQVMgcmVnaW9uYWxfcG9vbF9jb3VudCwKICAgICAgICAgICAgICAgICAgY291bnRpZihwcmVzZW50X2luX2NpdHkgPSB0cnVlKSBBUyBjaXR5X2NvdW50CiAgICAgICAgICAgICAgRlJPTSBtb2RlbC5yZWdpb25hbF9zcGVjaWVzX2ZpbHRlcmVkCiAgICAgICAgICAgICAgSk9JTiBtb2RlbC50YXhvbm9teSBVU0lORyAoc2NpZW50aWZpY19uYW1lKQogICAgICAgICAgICAgIFdIRVJFICMjQ09MX1NFTEVDVCMjID0gdHJ1ZQogICAgICAgICAgICAgIEdST1VQIEJZIHNjaWVudGlmaWNfbmFtZQogICAgICAgICAgKQogICAgICAgICAgU0VMRUNUIAogICAgICAgICAgICAgIHNwZWNpZXMuKiwKICAgICAgICAgICAgICBib2R5X21vcnBoc3BhY2UucGMxLAogICAgICAgICAgICAgIGJvZHlfbW9ycGhzcGFjZS5wYzIsCiAgICAgICAgICAgICAgYm9keV9tb3JwaHNwYWNlLnBjMywKICAgICAgICAgICAgICBib2R5X21vcnBoc3BhY2UucGM0LAogICAgICAgICAgICAgIHRyb3BoaWNfbmljaGUsCiAgICAgICAgICAgICAgZm9yYWdpbmdfbmljaGUsCiAgICAgICAgICAgICAgaXNfbm9jdHVyYWwsCiAgICAgICAgICAgICAgdHRjLmNsdXN0ZXIsCiAgICAgICAgICAgICAgdHRjLmNsdXN0ZXJfaGFsZiwKICAgICAgICAgICAgICB0dGMuY2x1c3Rlcl9kb3VibGUKICAgICAgICAgIEZST00gc3BlY2llcwogICAgICAgICAgSk9JTiBtb2RlbC50YXhvbm9teSBVU0lORyAoc2NpZW50aWZpY19uYW1lKQogICAgICAgICAgSk9JTiBkYXRhc2NpLnRheG9ub21pY190cmFpdF9jbHVzdGVycyB0dGMgVVNJTkcgKHNjaWVudGlmaWNfbmFtZSkiLCAnIyNDT0xfU0VMRUNUIyMnLCBjb2x1bW5fbmFtZSkKICAKICAgIHRiIDwtIGJxX3Byb2plY3RfcXVlcnkocHJvamVjdGlkLCBzcWwpCgogICAgZGF0YSA8LSBicV90YWJsZV9kb3dubG9hZCh0YikKICAgIHdyaXRlX2NzdihkYXRhLCBmaWxlbmFtZSkKICB9CiAgCiAgZGF0YSA8LSByZWFkX2NzdihmaWxlbmFtZSkKICAKICBkYXRhW2lzLm5hKGRhdGEkZm9yYWdpbmdfbmljaGUpLF0kZm9yYWdpbmdfbmljaGUgPC0gJ09tbml2b3JlJwogIAogIGRhdGEkZm9yYWdpbmdfbmljaGUgPSBhcy5mYWN0b3IoZGF0YSRmb3JhZ2luZ19uaWNoZSkKICBkYXRhJHRyb3BoaWNfbmljaGUgPSBhcy5mYWN0b3IoZGF0YSR0cm9waGljX25pY2hlKQogIGRhdGEkaXNfbm9jdHVyYWwgPSBhcy5mYWN0b3IoZGF0YSRpc19ub2N0dXJhbCkKICBkYXRhJHByZXNlbmNlX3JhdGlvID0gZGF0YSRjaXR5X2NvdW50IC8gZGF0YSRyZWdpb25hbF9wb29sX2NvdW50CiAgCiAgZGF0YQp9CmBgYAoKCmBgYHtyfQptZXJsaW5fZGF0YSA8LSBsb2FkX3NwZWNpZXNfZGF0YSgnbWVybGluJykKbWVybGluX2RhdGEKYGBgCgoKCmBgYHtyfQptZXJsaW5fcmVzdWx0IDwtIGRyZWRnZV9hbmRfc3Vic2V0KG1lcmxpbl9kYXRhKQptZXJsaW5fc3VtbSA8LSBtb2RlbF9zdW1tYXJ5KCdtZXJsaW4nLCAnc3BlY2llcycsIG1lcmxpbl9yZXN1bHQpCm1lcmxpbl9zdW1tCmBgYAoKYGBge3J9CmJpcmRsaWZlX2RhdGEgPC0gbG9hZF9zcGVjaWVzX2RhdGEoJ2JpcmRsaWZlJykKYmlyZGxpZmVfcmVzdWx0IDwtIGRyZWRnZV9hbmRfc3Vic2V0KGJpcmRsaWZlX2RhdGEpCmJpcmRsaWZlX3N1bW0gPC0gbW9kZWxfc3VtbWFyeSgnYmlyZGxpZmUnLCAnc3BlY2llcycsIGJpcmRsaWZlX3Jlc3VsdCkKYmlyZGxpZmVfc3VtbQpgYGAKCgpgYGB7cn0KYm90aF9kYXRhIDwtIGxvYWRfc3BlY2llc19kYXRhKCdib3RoJykKYm90aF9yZXN1bHQgPC0gZHJlZGdlX2FuZF9zdWJzZXQoYm90aF9kYXRhKQpib3RoX3N1bW0gPC0gbW9kZWxfc3VtbWFyeSgnYm90aCcsICdzcGVjaWVzJywgYm90aF9yZXN1bHQpCmJvdGhfc3VtbQpgYGAKCgpgYGB7cn0KZWl0aGVyX2RhdGEgPC0gbG9hZF9zcGVjaWVzX2RhdGEoJ2VpdGhlcicpCmVpdGhlcl9yZXN1bHQgPC0gZHJlZGdlX2FuZF9zdWJzZXQoZWl0aGVyX2RhdGEpCmVpdGhlcl9zdW1tIDwtIG1vZGVsX3N1bW1hcnkoJ2VpdGhlcicsICdzcGVjaWVzJywgZWl0aGVyX3Jlc3VsdCkKZWl0aGVyX3N1bW0KYGBgCgotLS0tLS0tLS0tLS0tLS0tLS0tLS0tCkZ1bGwgcmVzdWx0Ci0tLS0tLS0tLS0tLS0tLS0tLS0tLS0KYGBge3J9CmFsbF9zcGVjaWVzX3Jlc3VsdHMgPC0gZnVsbF9qb2luKGZ1bGxfam9pbihtZXJsaW5fc3VtbSwgYmlyZGxpZmVfc3VtbSksIGZ1bGxfam9pbihib3RoX3N1bW0sIGVpdGhlcl9zdW1tKSkKd3JpdGVfY3N2KGFsbF9zcGVjaWVzX3Jlc3VsdHMsICJzcGVjaWVzX3Jlc3VsdC5jc3YiKQpgYGAKCmBgYHtyfQptZXJsaW5fZGF0YSRzb3VyY2UgPC0gJ01lcmxpbicKYmlyZGxpZmVfZGF0YSRzb3VyY2UgPC0gJ0JpcmRsaWZlJwpib3RoX2RhdGEkc291cmNlIDwtICdCb3RoJwplaXRoZXJfZGF0YSRzb3VyY2UgPC0gJ0VpdGhlcicKCnBsb3RfZGF0YSA8LSByYmluZChtZXJsaW5fZGF0YSwgYmlyZGxpZmVfZGF0YSwgYm90aF9kYXRhLCBlaXRoZXJfZGF0YSkKCmdncGxvdChwbG90X2RhdGEsIGFlcyh4ID0gcGMxLCB5ID0gcGMyLCBhbHBoYSA9IHByZXNlbmNlX3JhdGlvLCBjb2xvciA9IGFzLmZhY3Rvcihmb3JhZ2luZ19uaWNoZSkpKSArIGdlb21fcG9pbnQoKSArIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikgKyBmYWNldF93cmFwKH4gc291cmNlKQpgYGAKYGBge3J9CmdncGxvdChwbG90X2RhdGEsIGFlcyh4ID0gcGMzLCB5ID0gcGM0LCBhbHBoYSA9IHByZXNlbmNlX3JhdGlvLCBjb2xvciA9IGFzLmZhY3Rvcihmb3JhZ2luZ19uaWNoZSkpKSArIGdlb21fcG9pbnQoKSArIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikgKyBmYWNldF93cmFwKH4gc291cmNlKQpgYGAKYGBge3J9CmxpYnJhcnkoZ2dwdWJyKQpgYGAKCmBgYHtyfQptZXJsaW5fZGF0YSRuaWNoZSA8LSBwYXN0ZShtZXJsaW5fZGF0YSR0cm9waGljX25pY2hlLCAnICcsIG1lcmxpbl9kYXRhJGZvcmFnaW5nX25pY2hlLCAnICcsIG1lcmxpbl9kYXRhJGlzX25vY3R1cmFsLCAnICcsIG1lcmxpbl9kYXRhJGNsdXN0ZXIpCmBgYAoKYGBge3J9Cm15cGxvdCA8LSBnZ2FycmFuZ2UoCiAgbmNvbCA9IDIsCiAgbnJvdyA9IDIsCiAgZ2dwbG90KG1lcmxpbl9kYXRhLCBhZXMoeCA9IHBjMSwgeSA9IHBjMywgYWxwaGEgPSBwcmVzZW5jZV9yYXRpbywgY29sb3VyID0gbmljaGUpKSArIGdlb21fcG9pbnQoKSArIHRoZW1lX2J3KCkgKyAgdGhlbWUobGVnZW5kLnBvc2l0aW9uPSJub25lIiwgYXhpcy50aXRsZS54PWVsZW1lbnRfYmxhbmsoKSkgKyB5bGFiKCJTVCArIFBCIC0+IExUICsgU0IiKSwKICBnZ3Bsb3QobWVybGluX2RhdGEsIGFlcyh4ID0gcGMyLCB5ID0gcGMzLCBhbHBoYSA9IHByZXNlbmNlX3JhdGlvLCBjb2xvdXIgPSBuaWNoZSkpICsgZ2VvbV9wb2ludCgpICsgdGhlbWVfYncoKSArIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbj0ibm9uZSIsIGF4aXMudGl0bGUueD1lbGVtZW50X2JsYW5rKCksIGF4aXMudGl0bGUueT1lbGVtZW50X2JsYW5rKCkpLAogIGdncGxvdChtZXJsaW5fZGF0YSwgYWVzKHggPSBwYzEsIHkgPSBwYzQsIGFscGhhID0gcHJlc2VuY2VfcmF0aW8sIGNvbG91ciA9IG5pY2hlKSkgKyBnZW9tX3BvaW50KCkgKyB0aGVtZV9idygpICsgdGhlbWUobGVnZW5kLnBvc2l0aW9uPSJub25lIikgKyB5bGFiKCJMVCArIFBCIC0+IFNUICsgU0IiKSArIHhsYWIoIkJvZHkgU2l6ZSIpLAogIGdncGxvdChtZXJsaW5fZGF0YSwgYWVzKHggPSBwYzIsIHkgPSBwYzQsIGFscGhhID0gcHJlc2VuY2VfcmF0aW8sIGNvbG91ciA9IG5pY2hlKSkgKyBnZW9tX3BvaW50KCkgKyB0aGVtZV9idygpICsgdGhlbWUobGVnZW5kLnBvc2l0aW9uPSJub25lIiwgYXhpcy50aXRsZS55PWVsZW1lbnRfYmxhbmsoKSkgKyB4bGFiKCJCZWFrIFNpemUiKQopCm15cGxvdApgYGAKYGBge3J9CmpwZWcod2lkdGggPSAxMDI0LCBoZWlnaHQgPSAxMDI0KQpteXBsb3QKZGV2Lm9mZigpCmBgYA==